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We describe recent progress in the numerical study of the structure of rapidly 
rotating superfluid neutron star models in full general relativity. The superfluid 
neutron star is described by a model of two interpenetrating and interacting 
fluids, one representing the superfluid neutrons and the second consisting of 
the remaining charged particles (protons, electrons, muons). We consider gen¬ 
eral stationary configurations where the two fluids can have different rotation 
rates around a common rotation axis. The previously discovered existence of 
configurations with one fluid in a prolate shape is confirmed. 


1 Introduction 

Most current studies of neutron star oscillations are based on a simple perfect 
fluid model for a neutron star, which neglects (besides other things like the magnetic 
field, crust etc.) the crucial importance of superfluidity. The presence of substantial 
amounts of superfluid matter in neutron stars is backed by a number of theoreti¬ 
cal calculations of the state of matter at these extreme densities (e.g. see |LjJ, j^]), 
and by the qualitative success of superfluid models to accommodate observed fea¬ 
tures of glitches and their relaxation (e.g. see U and references therein). For both 
the study of oscillation modes of superfluid neutron stars, and the modelling of the 
particular glitch instability, it is important to know the (quasi-) stationary initial 
state of the unperturbed neutron star. The principal difference with ordinary-fluid 
models consists of the fact that in a superfluid model there are two fluids (i.e. the 
superfluid neutrons and all the rest), that can flow independently. As a result of the 
slow spindown process of the neutron star, the neutron superfluid will lag behind the 
charged constituents, due to its complete absence of viscosity. The generic configu¬ 
ration therefore corresponds to a (quasi-) stationary two-fluid model, with the two 
(strongly coupled!) fluids rotating around the same axis, but with different rotation 
rates. This model has been studied before in both the Newtonian (i!) and the 
relativistic framework ||, using the slow-rotating approximation and a simple “poly¬ 
tropic” two-fluid equation of state. Here we describe a fully relativistic numerical 
approach to solve the problem for arbitrary rotation rates. 


1 


2 


R. Prix, J. Novak and G.L. Comer 


2 The two-fluid model 


We briefly present the notation and formalism used to describe the two-fluid 
model for superfluid neutron stars, a more detailed discussion can be found in fj], |j . 
We consider two fluids, representing the neutrons (denoted by the label n) and the 
remaining charged components (denoted by c) respectively. Superfluidity allows the 
neutrons to flow independently of the remaining constituents, without any direct 
frictional force between the neutrons and the second fluid. Another implication 
of superfluidity is that the rotating neutron fluid will be threaded by a lattice of 
microscopic quantized vortices, but their effect on the macroscopic scale is neglected 
here, so we treat the neutrons simply as an inviscid fluid. The kinematics of the two 
fluids is described by the respective rest-frame particle number densities n n and n c , 
together with the four-velocities tt„ and u c . The corresponding particle four-currents 
are therefore n„ = n n u n and n c = n c u c , and we assume strict conservation of both 
fluids (i.e. we exclude the possibility of “transfusion” || via /3-reactions between 
protons and neutrons), which means 


V • n n = 0 , and V • n c = 0 . 


( 1 ) 


The dynamics is governed by the total energy density £(n B , n c , A 2 ), i.e. “equation of 
state”, where A represents the relative velocity between the two fluids, and is defined 
as 


A = 1 - 


n n n c 


( 2 ) 


n n - n 

corresponding to a gamma factor T A of the relative motion, which is 

r A = ~u n ■ w c = (l - a 2 ) _1/2 . (3) 

The energy scalar £ determines the first law of thermodynamics, namely 

d£ = p n dn n + p c dn c + a dA 2 , (4) 


which defines the chemical potentials p n and /A of the two fluids, and the “entrain¬ 
ment function” a, which measures the dependence of the energy density on the rela¬ 
tive velocity. The equations of motion for the two interacting fluids are derived from 
a “convective” variational principle || based on the hydrodynamic Lagrangian 
£h = —£■ This defines the canonical momenta p n and p c via the total differential 
5 Ah = P n ■ 5n n + p c ■ 5n c . These canonical momenta are found explicitly as 


2a 


p — /i u n T ^ (u c r^w n ), 

a 

r r 2« . 

p = P U c -\ - -T (tin - r A tt c ) , 


™ c Ia 


( 5 ) 

( 6 ) 
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which shows that the effect of the entrainment a is to make the momenta deviate 
from their respective four-velocity directions. This effect has first been discussed (in 
different terms) by Andreev and Bashkin in the context of mixtures of superfluid 4 He 
and 3 He 0> and is also equivalent to a description in terms of “effective masses” 
[||. In the absence of mutual friction and pinning forces between the two fluids, the 
equations of motion have a very simple form, despite the fact that two fluids are 
coupled via the equation of state, namely 

u n ■ dp n = 0 , and u c ■ dp c = 0 , (7) 

where d denotes the exterior derivative, i.e. ( dp ) a 0 = V a P 0 — Vpp a . The corre¬ 
sponding two-fluid energy-momentum tensor T is obtained from the variational prin¬ 
ciple as 

T = n a (g> p n + n c 0 p c + Tg , (8) 

where g is the metric tensor, and 'k is the generalized pressure, defined as 

\k + £ = n n p n + n c g° . (9) 


3 Symmetries and integrals of motion 

We consider spacetimes that are stationary, axisymmetric and asymptotically 
flat. Stationarity and axisymmetry are associated with two (commuting) Killing vec¬ 
tor fields t and ip respectively. We choose adapted polar-type coordinates {t, r, d, p}, 
which are such that t = d/dt and ip = 8 /dp. We assume the two fluids are in purely 
circular motion around the 2 -axis (i.e. we exclude convective meridional currents), 
so we write 

u a = 7 n {t + S! n <p) , and u c = 7 c (t + fi c tp) , ( 10 ) 

where fl n and H c are the respective rotation rates of the two fluids. The most general 
metric (in quasi-isotropic maximal slicing gauge) [pd]] can be written as 

ds 2 = -(N 2 - N V N V ) dt 2 - 2N v dpdt + A 2 (dr 2 + r 2 d0 2 ) + B 2 r 2 sin 2 6 dp 2 , (11) 

where N, NA and B are functions of (r, 0), and N v = g w N v . With the additional 
assumption of uniform rotations (Vf2 = 0), the equations of motion © simplify 
greatly and can be shown to reduce to two first integrals of motion ||, namely 

p n = 7 „ const n , and p c = y c const c . ( 12 ) 

4 Results and discussion 

Under the assumptions of the previous section, the Einstein equations reduce 
to four coupled elliptic equations for the four unknown functions N, N' p , A and 
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B, with the source terms given by the energy-momentum tensor (||). These can 
be solved very accurately using an iteration scheme involving the first integrals of 
motion ©• The numerical scheme is based on a multi-domain pseudo-spectral 
method, which is described in more detail in [|l2j, and which we have appropriately 
extended to the two-fluid model considered here. In principle the problem can be 
solved for any equation of state (EOS), but currently only a polytropic toy-model 
EOS is implemented. This, however, is sufficient to show the qualitative features of 
a superfluid neutron star models as opposed to a single fluid model. The EOS used 
here is a “generalized polytrope” of the form 

£ = ( m n n n + m c n c )c 2 + ^K n n? + ^c^c 2 + K nc ril 3 n2 4 + f3n2 5 n2 e A 2 . (13) 

We note that the two fluids are coupled via a “symmetry energy”-type term propor¬ 
tional to K nc and an “entrainment” term proportional to /3A 2 . 

The numerical two-fluid code is currently in the test phase, internal consistency 
checks like the relativistic viriel (cf. (TTJ) are satisfied to about 10 -5 at a convergence 
of 10 -6 . However, further tests are necessary, in particular comparisons with previous 
results in the slow-rotation approximation, both relativistic Q and Newtonian Q, ||. 
The latter, in particular, allowed an analytic solution for a subclass of the EOS 
( p~3| ) . Furthermore, the existence of prolate solutions for the slower rotating fluid was 
found analytically in Q, and we have been able to reproduce these solutions with 
our numerical code (see Fig. []]). In Figure |l| we see the influence of the interaction 
(characterized by K nc and (3) between the two fluids on the stationary configuration 
of a neutron stars model with roughly M ss 1.4 Mq and R ss 11 km. The EOS 
parameters used are 7i = 72 = 2, and 73 = 74 = 75 = 76 = 1, and K n = k c = 0.06. 
In this configuration, one fluid is rotating at fl n = 2000 7r s -1 , while the second 
fluid rotates much slower at only fl c = 200 7r s -1 . The left-hand graph shows the 
solution for two fluids which are only gravitationally coupled (i.e. K nc = j3 = 0), 
while in the right-hand graph there is additional coupling via the “symmetry”-term 
Knc and entrainment f3. The effect is to push the faster fluid further out and to 
“squeeze” the slower fluid, in this case to the point of actually making it prolate. 
Both the uncoupled and the coupled configuration are characterized by the same 
central chemical potentials /i n (0) = /U c (0) = 0.3 me 2 , which is why their resulting 
total masses and radii differ by about 10 %. 

Apart from further quantitative tests, more work is necessary in order to allow 
for more realistic equations of state, and to include the possibility of vortex-mediated 
pinning forces between the two fluids. 
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Figure 1: Two-fluid configurations with O n = 2000 n s _1 and fl c = 200 n s _1 . The 
thick lines represent the surfaces of the two fluids respectively, while the dashed lines 
are the isosurfaces of “gravitational potential” In N. In the left-hand configuration, 
the two fluids are locally non-interacting, i.e. k„ c = 0 and (3 = 0, while in the 
right-hand configuration we have chosen n nc = 0.02 and (3 = 0.12. 
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